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ABSTRACT 


This work represents the results of one phase of research conducted 
for the JPL Solar Electric Propulsion (SEP) Navigation Software System 
development program. It deals only with the problem of designing the flight 
quality trajectory program, which is a major subset of the entire navigation 
software system. 

In this phase of research (breadboard development phase), attempts 
were made to assess the SEP trajectory software functional requirements, 
to investigate the program design method satisfying these requirements, to 
identify the primary anticipated problem areas, and to provide solutions to 
these problem areas. These efforts culminated in the development of a 
compact breadboard program, "LOWTRAJ. " A functional description and 
the mathematical formulation of the program are presented. 

The results of tests performed using LOWTRAJ indicate that the 
primary requirements of the flight quality trajectory program can be met 
with this type of design. Future extensions of the program, further refined 
to support flight operation, should be straightforward. 
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I. INTRODUCTION 


In the past decade, various trajectory programs have been developed 
to study the applications of solar electric propulsion (SEP) in interplane- 
tary explorations. These programs were of a first-generation type, 
oriented primarily to support the assessment of SEP mission feasibility, 
payload performance, and reliability. In this early conceptual stage of SEP 
technology, the hardware performance characteristics were not well defined. 
Simplified approximations and idealizations were used to represent low- 
thrust propulsion, which made the mathematical formulation and implementa- 
tion of optimal trajectory search feasible. This resulted in a large quantity 
of preliminary information stressing the upper limit of SEP spacecraft per- 
formance capabilities. It appears that with currently available technology a 
SEP spacecraft can be an attractive candidate for some selected missions. 

To prepare for an actual SEP flight, second-generation SEP trajectory 
software is being considered. It is intended to be a low-cost program, 
serving primarily as a test tool (breadboard) for the development of a third- 
generation flight program. It is also intended to become the main core of a 
larger breadboard SEP Navigation Software System (SEPNSS), which will 
include the trajectory, orbit determination (OD), and guidance control 
software. 

For the support of a SEP flight, particularly in performing reliable 
navigation and guidance functions, it is critical to have a very accurate 
"theoretical prediction" of the spacecraft state. An accurate simulation of 
SEP spacecraft controls and their net effect is required to achieve this goal. 
Optimum payload performance to the last few percent, the characteristic of 
the first-generation software, is now considered to be of secondary 
importance. To design a practical mission, optimality must be violated 
anyway to accommodate specific restrictions required in hardware design, 
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and to satisfy other mission constraints imposed for the sake of mission 
reliability, scientific experiments, and communications. The development 
of elegant mathematical approaches using the calculus of variations to 
generate suboptimal controls satisfying many and varied types of constraints 
appear to be rather difficult, and the potential performance margin to be 
gained is expected to be slight. Therefore, a direct parametric trajectory 
search procedure that may or may not explicitly optimize the payload is 
proposed. In addition to satisfying the requirements of establishing a valid 
reference trajectory, the trajectory program is expected to interact 
intimately with the OD and guidance programs. The use of the direct 
parametric search procedure is most natural, and in harmony with the needs 
of these user programs. 

In Section II, a general description of a trajectory program is given. 
Details of specific problems and solutions required for the low-thrust 
application will be discussed in Section III. Section IV will be devoted 
entirely to the mathematical formulations and solutions of these specific 
problems. 


II. GENERAL TRAJECTORY SOFTWARE FUNCTIONS 
AND REQUIREMENTS 

A simplified diagram of a SEPNSS is shown in Fig. 1. This illustrates 
typical (ballistic or SEP) OD and guidance software functions, and their 
relationships to one another. The subprograms PATH, VARY, and SEARCH 
are the major structural constituents of the LOWTRAJ. As one can infer 
from the diagram, the design philosophy of the OD and guidance programs 
must strongly influence the formulation of PATH and VARY. Since it is 
impractical to expect synchronized progress in the development of the OD, 
guidance, and trajectory programs, the major linkages of PATH- VARY with 
the OD or guidance programs will not be attempted until later in the SEPNSS 
development phase. 

Descriptions of the general functions and requirements for each of the 
three major subprograms of LOWTRAJ are given in the following. Require- 
ments unique to the low-thrust applications, demanding special study, are 
emphasized. 
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A. 


Trajectory Simulation Program (PATH) 


The function of this program is to theoretically predict the spacecraft 
path by numerical integration of the equations of motion. It must accurately 
account for and model all forces acting on the spacecraft, particularly the 
low- thrust for ce , and must integrate , maintaining high numerical accuracy. 
This is an open loop path predictor that does not perform a targeting 
function. 

B. Variational Equations Integrator (VARY) 

This program generates variations of the spacecraft trajectory as 
induced by small perturbations of the trajectory parameters {P} (i.e. , 
9X(t)/3P). This is accomplished by the integration of the variational 
equations. These are the basic data required by the user programs, namely 
trajectory search, OD, and the guidance programs. The major objective 
in the construction of this subprogram is to identify the important param- 
eters with respect to which partial derivatives are required by the search 
program. At the same time the potential needs of the OD and guidance 
programs must be considered. Then the derivation and integration of the 
necessary variational equations must be performed, and the efficient 
transfer of this information to the user programs must be executed. 

C. Trajectory Search Program (SEARCH) 

This program performs the deterministic targeting function. It drives 
PATH and VARY, and iteratively searches for the trajectory shaping 
parameters that will satisfy the required boundary conditions and various 
mission constraints. A large number of search parameters (perhaps as 
many as 100 or more) can be involved in this process. This program must 
handle searches with many degrees of freedom, preferably with a capability 
to optimize the payload. It must have the ability to handle many constraints 
of different types, and it must attain convergence with high reliability and 
speed. 
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III. SPECIFIC PROBLEMS AND PROPOSED SOLUTIONS 


A. Accurate Low- Thrust Modeling in PATH 

An accurate modeling of the low-thrust force depends on a good 
understanding of the performance characteristics of the thrust-producing 
hardware subsystems, the policies of their operation, and the many modes 
in which they may be operated. Mathematical simulations of these control 
processes by means of simple equations and parameterization leads to 
correct representation of the net force exerted on the spacecraft. The key 
to successful modeling lies in the manner in which these parameters are 
introduced. The versatility of the program, its compatibility to the search, 
OD, and guidance programs are established at this point. 

In the following, a brief discussion of the three major thrust-producing 
hardware systems, and the assumptions being made on their operating 
policies and modes will be presented. The resulting model and the meaning 
of some important parameters will be analyzed. 

1. Simulation of the power subsystem function. Power subsystem 
components pertinent to low-thrust control include the solar array and 
switching and control function subsystems for power management and dis- 
tribution. Theoretically, solar array maximum output power minus house- 
keeping power should be available to the thrust subsystem. However, a 
closed loop maximum power operating policy for the thrust subsystem may 
not be acceptable. Unpredictable fluctuations in solar array output power, 
compounded by the already noisy thruster performance at known operating 
levels, would make the overall low-thrust noise level too large to be 
tolerable for accurate navigation and guidance. Therefore, regardless of 
the maximum power available for propulsion from the solar array, the 
power input to the thrust subsystem will be programmed to be less than 
the maximum point. The devices for the solar array power regulation 
have as yet to be specified by the power subsystem specialists. Still, the 
regulation generally consists of triggers to command the level changes, and 
automatic maintenance of the set level between the triggers. The on-board 
command and control timer triggers these changes, and the interval between 
the triggers is estimated to be on the order of 1 to 10 days depending on the 
power profile. Basically, the power regulation can be achieved in two ways. 


4 


JPL Technical Memorandum 33-648 


First, it is assumed that the regulation device automatically controls both 
the current and the voltage output of the solar array at a fixed value between 
regulation switching times (for convenience call it a "power stage 11 ). This 
implies a piecewise constant power operation. Second, if the regulation 
device controls only the current (or the voltage) and lets the voltage (or the 
current) operate at the natural solar array output, the program must model 
the solar panel cur rent- voltage (I-E) characteristic curves as a function of 
solar distance. This is not only more complex in mathematical formulation, 
but it also adds greater uncertainty in the thrust magnitude, because the 
accuracy of the given (I-E) curves as a function of solar distance is doubtful. 
If there are any uncertainties, these must be fed to a navigation program 
and their impact measured. The present trajectory software package 
includes only the first policy and models power input to the thrust subsystem 
as being piecewise constant. 

The lower bound of solar array maximum output power during t . to 

t[ i (i-th power stage) is estimated based on the spacecraft state (r, *r) at 

time t. and the conventional power curve formula. This power minus the 

housekeeping power is the available input power to the thrust subsystem at 

the i-th stage [P (i)]. Actual input power to the thrust subsystem will be 
a 

[v.P (i)], which comprises the basis for panel power regulation. Here v. 

1 a i 

is a power utilization factor for the i-th power stage, = 1 represents full 
utilization, = 0 represents a coast period, and 0 < v\ < 1 represents 
partial utilization. Even though a bang-bang-type control, where = 0, or 
1, is the optimal control policy, one may want to design the nominal path 
with v\'s slightly less than 1 to provide some guidance reserve. In addition 
v\ T s can be used to simulate expected or unpredictable solar array degrada- 
tions caused by solar flares or meteorite impacts. It may even be assigned 
a value larger than 1 to simulate conditions where the actual output of the 
panel in space indeed exceeds the theoretical prediction. 

2. Simulation of the thrust subsystem function. The two major 
thrust subsystem components considered are pbwer conditioner (PC) units 
and thruster (THR) units. A combined operation of one PC unit and a 
thruster constitutes a thrust unit. 
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The characteristic numbers of the thrust subsystem are the 

maximum and minimum power ratings of a thrust unit, the number of thrust 

units available, the efficiency of power conversion, and the specific impulse 

(I ). The efficiencies of the PC and THR depend on the respective input 
sp 

power levels. The specific impulse depends on the operating power 
(throttling level) of the THR. Currently available data on these (Ref. 1) can 
be represented adequately by quadratic functions of the input power to the 
thrust unit. 

The operating policy of the thrust subsystem assumes that the total 

input power to the thrust subsystem [V.P (i)] is nominally distributed 

equally among the minimum number (N m ^ n (i)) of thrust units required to 

match the power. This is a maximum efficiency policy for a given power 

level. If the number of available thrust units are less than the N . (i), then 

min 

all units will be operated at maximum level. As an option one may specify 
the minimum number of units to be operated. This could, for example, be 
the case if one wants to maintain at least two thrusters operating, so that 
three-axis attitude control can be maintained using low thrust. 

3. Simulation of the thrust vector control system function. Thrust 
vector control is assumed to be implemented with the aid of sun and star 
sensors for attitude reference. Gross reorientation of the thrust vector is 
achieved by gimballing the sensors by a desired amount, thus offsetting the 
tracking, then applying torque to the spacecraft to reacquire the stars. 

This maneuver will be performed at specified intervals (call it angular 
stages), which may range from a few to hundreds of days depending on the 
mission. Between gross reorientations, the autonomous attitude control 
command system maintains the thrust vector cone and clock angles within a 
specified tolerance band. 

B. Setup for Variational Equations in VARY 

The dimensions of the partial derivatives to be given by VARY 
depend on the needs of the user programs. SEARCH requires partial 
derivatives with respect to injection state, solar array power at 1 AU 
(Pq), thrust vector directions (two angles) for all angular stages, power 
utilization factors ^ for all power stages, the arrival time, and the arrival 
velocity bias of the spacecraft with respect to the target. 
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Further expansion of VARY to meet the needs of OD and guidance 
functions is as yet to be determined. 

C. Design of Program SEARCH 

1, Input requirements. As discussed previously, this program is 
of a second-generation type emphasizing the accuracy in the trajectory 
through accurate hardware function modeling. The capability to search for 
optimal gross mission parameters such as Pq, 1^, launch dates, and flight 
time is considered to be outside the scope of this program. Nevertheless 
optimal or at least near-optimal mission profile is preferred. Therefore, 
the first guess of the general thrust profile is obtained from the first- 
generation trajectory optimizing programs, such as CHEBYTOP or EPITOP 
(Refs. 2, 3). This is mandatory not only for the sake of performance, but 
also for easier convergence. This program does not have the capability to 
self- start, nor is it meant to generate the "ballistic conic path" equivalent 
of a low-thrust trajectory. Once the crude profile is given, it will readjust 
all the free search variables by means of a modified Newton Raphson 
method to satisfy the required boundary conditions, one of which may include 
the final mass with a given tolerance of, say, plus or minus 5 to 10 kg. It 
can, as an option, perform a limited optimal search for the maximum final 
mass. 


2. Versatility. The flexibility of the program is specifically geared 
for the needs of flight project analysis and design. It must accurately 
simulate the controls of specific form required for mission implementation. 
It must also be able to simulate various types of control malfunctions so that 
impacts of these uncertain hardware functions can be analyzed and a reliable 
mission designed. 

These goals can be achieved if one allows, by option, all the low- 
thrust control parameters to be included in the search or to be fixed. The 
only parameters that cannot be given this freedom are the stage times, both 
for power stages or angle stages. However, the user will have the freedom 
to assign almost any stage pattern as long as it does not exceed the desig- 
nated dimension of the stages in the program, which currently is 200 for 
power stages and 50 for angle stages. 
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Some examples of the desirable features are: 

(1) The user can design the trajectory with thrust on-off time 
specified. In first generation trajectory software, this was 
not possible because the optimality condition controlled and 
internally generated the switching time. Coast phases could 
not be arbitrarily specified. This capability is crucial; it will 
be needed to satisfy the science experiment requirements, navi- 
gational needs, and for reliable mission design. 

(2) Since all the power utilization parameters can be searched or 
fixed, any throttling levels may be commanded. 

(3) Thrust vector and/or spacecraft attitude can be constrained 
for any desired period of time. Such constraints are imposed 
usually by the limited structural flexibility of the spacecraft, 
thrust-subsystem thermal control requirements, the science 
experiments or the communications requirements. 

(4) Thruster arcing or failures can be simulated and updates of 
trajectory can be made. 

(5) Solar panel degradations, minor or major, can be simulated 
and their impacts can be measured. 

3. Organization of search variables. 

(1) Independent variables: To attain maximum flexibility, a large 

dimension in independent-variable space is introduced. This 
high degree of freedom consists mainly of thrust angles modeled 
in multistage fashion. Since SEP is a power-limited propulsion 
device, these angles are the main source of trajectory shaping 
capability, particularly after spacecraft initial injection. 

Other important degrees of freedom that can be used for control 
are the thrust duration and the time of encounter. The search 
on thrust duration is performed on allowing only quantum 
jumps, that is 0 to 1, or 1 to 0. 

The degrees of freedom of 0, 2, and 3 can be assigned, by 
option, to departure and arrival velocity biases (Vg). 
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These correspond to cases where it does not search on Vg 
(rendezvous), Vg is constrained in magnitude but has 2 degrees 
of freedom in the choice of direction (flyby with given relative 
speed), and, thirdly, Vg is unconstrained. 

Other variables such as injection mass, injection time, and 
P Q are included to meet the primary needs of flight quality 
mission design. 

(Z) Dependent variables: The program is organized in a manner 

such that dependent variables always include the spacecraft final 
position and velocity minus position and velocity biases. These 
are always searched to coincide with the state of the target body. 
Analytic ephemerides of the major planets, asteroids, and 
comets are internally linked to the program. 

An additional dependent variable included is the final mass of 
the spacecraft. This is included so that limited optimal control 
of the final mass can be accomplished. 

4. Search procedure . As it has been stressed in Section III, the 
many degrees of freedom of search are due to the many angle variables. 
Since these are characteristically the same controls appearing consecutively 
and progressively in stages, it is likely that high correlations exist among 
the partial derivatives [BY/B(^j, Pj)] (where Y is the vector of dependent 
variables, and Q'j and pj are the thrust cone and clock angles for stages 
1=1,- • • , N, etc. , to be searched). It is unwise to perform standard 
iterative linear analysis of the form MaX = AY without fully analyzing the 
singularity of matrix M. The search algorithm makes extensive utilization 
of the information obtained in performing the "Singular- Value Decomposition 
Analysis" (Ref. 4 ) of matrix M. Methods of obtaining solutions within the 
framework of linear algebra is discussed in greater detail in Section IV. 


JPL Technical Memorandum 33-648 



IV. MATHEMATICAL FORMULATION OF LOWTRAJ 


f 


A. Spacecraft Trajectory 

Generally, the sources of acceleration of the spacecraft relative to the 
center of integration to be considered include the following perturbations: 

(1) The Newtonian point-mass acceleration relative to the center of 
integration. 

(2) The acceleration due to low thrust. 

(3) The acceleration due to chemical motor burns. 

(4) The acceleration due to solar radiation pressure. 

(5) The accelerations due to other smaller order gravitational 
interactions, including N-body effects, planet oblateness effects, 
mascon effects, and relativistic effects. 

(6) The accelerations due to small perturbations originating in the 
spacecraft, attitude controls (especially the low-thrust type), 
and due to gas leaks. Nonavailability of solar power for low 
thrust during solar occultation must also be included. 

However, due to the experimental nature of the program and to maintain low 
cost and efficiency of program operation, numerical integration is performed 
in single precision. Inclusion of perturbations (4) to (6) in the single pre- 
cision algorithm is not meaningful, thus they are excluded in LOWTRAJ but 
will be required in the flight quality program. In the current scheme of the 
trajectory search program, the inclusion of these small forces is not 
expected to influence the basic algorithms. 
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B. 


Glossary of Notations 


Unless otherwise stated, the following notations will be used 
consistently without explanation. 

F low-thrust acceleration 

f low-thrust magnitude 

G gravitational acceleration 

H mass flow rate 

m spacecraft mass at time t 

mg. spacecraft mass at the beginning of i-th power stage 
7 spacecraft position vector at time t 

7^ spacecraft position vector at the beginning of i-th power stage 

7 unit vector of specified reference star position 
h initial time of power stage i 

Vq^ spacecraft velocity vector at the beginning of i-th power stage 
v spacecraft velocity vector at time t 

X spacecraft state vector at time t where X = (r , v, m) 

Xq^ spacecraft state vector at the beginning of i-th power stage 

a, (3 thrust cone and clock angles with respect to sun and a star 
At. i-th power stage interval 

Atj I-th angular stage interval 

(i gravitational constant of the sun 

v power utilization factor 

I low-thrust unit vector 
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Subscripts 

i power- stage 

I angular-stage 

C. Derivation of Equations of Motion 

1. Mathematical force model . Mathematical expressions of the 
spacecraft accelerations and mass flow rate are given as follows. 

a. Gravitational acceleration 

5(7) = (i) 

r 

b. Low-thrust acceleration 


F Ii (r 0i’ v 0i' °r ^1* v i ) “ 


f i ^Oi* V v i ) % {a r ?) 


m 


( 2 ) 


where L (?q*> v^, i^) is the thrust magnitude for the i-th power stage, and 

| (OTj, (3 , ?) is the thrust unit vector during I-th angle stage. 

Note: v. is constant during the i-th power stage and 

(3j are constants during the I-th angle stage. 

c. Mass flow rate 


H i " f i ^ r 0i’ v 0i* v i^ c 


where c is the thruster exhaust velocity. 
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Derivation of f. fr^, v^., v.) 
Let 


r 

P 


r 0i + 2 1 


r 0i ‘ V 0i 
r 0i ’ V 0i 


+ 1 


r 0i ' v 0i 


l r 0i' 


At. 

l 


(3) 


where r is the estimated upper bound of the spacecraft- sun distance, and 
P 

this is used to estimate the lower bound of the solar-panel maximum output 
power during At.. To maintain the validity of such a power estimate, AL 
must not be too large. Then, 


P 


max 


i - 1 


a.r 


i P 


- (i+3)/2 


(4) 


where P max = estimated lower bound of panel maximum output power during 
the i-th power stage, a^ = solar panel maximum power curve coefficients 
and Pq = solar panel output power at 1 AU. Let 

**a ~ Vax ~ **h 

where p^ = available input power to the thrust unit and p^ = housekeeping 
power. Let 


r v ‘ 

N . = - 1 + 1 (integer operation) (5) 

min p 

r l 


P v - 

N = — 1 (integer operation) (6) 

max p 

r 2 
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Then 


If N . > N then N = N 

min a o a 

N = N . or / If N . < N , then N = N 

o min \ min d o d 

If N > N then N = N 


max 


max 


where 


N . 
mm 


N 

max 



minimum number of thrust units required to be in operation 

maximum number of thrust units one may be operating 
without throttling below the minimum rated power 
maximum power rating of a thrust unit 

minimum power rating of a thrust unit 

number of thrusters actually operated 

available number of thrust units 

desired lower limit of thrust units operating 


Let 


P v - 
r a i 


op ' N 


C = C 0 + C 1 P op + C 2 Pop 


ti n 0 + P Q p + ^2 p op 


(7) 


( 8 ) 
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where 


P Q p = the operating level of each thrust unit 

c = thruster exhaust velocity (I^g) 

H = thrust unit power conversion efficiency 
Cq, Cj, c^ = polynomial coefficients used to express exhaust velocity 
as a function of operating power 

V ^2 = efficiency coefficients 


Then 


f. 

i 


2ri v. p 

1 ra 
c 


(9) 


Spacecraft mass flow rate 


H. 

l 


2n v.p 
1 la 


( 10 ) 


where fL = mass flow rate during i-th power stage. 


Derivation of £( a j> Pj, 7). Let 


k' = j‘ = k' X s, i' = k' X j' 


( 11 ) 


where s is the reference star unit vector and i', j', k' are unit vectors of 
sun- star reference frame. Let 


= (sin cos |3j, sin sin pj, cos a^) (12) 
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where £' is a thrust unit vector in the sun- star reference frame and a^, (3^. 
are the thrust cone and clock angles. Then, 

I = T(r’)f' (13) 


where T(r) is a coordinate transformation matrix with the following 
components: 



(14) 


where subscripts x, y, and z denote the x, y, and z components of the 
vector. 

Equations of motion. To maintain the symmetry in the expression, the 
equations of motion are expressed in seven first order differential equations. 


“ 


• ■ 



— 

r 


V 




V 


G + F t . 
Ii 

rh 


H. 



l 


Since the first derivatives X are discontinuous at the bounds of every power 
stage and angular stage (which are designed to coincide with one of the 
power-stage times), numerical integrations are performed piecewise by 
power stage increments with a restart procedure for each discontinuity. 
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D. Derivation of Variational Equations 


Partial derivatives of the spacecraft state with respect to search 
parameters are obtained by numerical integrations of the variational 
equations. In functional form, equations of motion are given by 

X = jf(x, Q) (16) 


where Q is the parameter set with respect to which partial derivatives are 

required. For search purposes, the parameter set Q includes (Xqj, 1J), 

where cf = (a , P v . , p ) for all i and I desired. 

1 1,10 


The variational equations are of the following general form: 





a^ix, Q) ax , a^r(x, Q) 


ax 


9Q 3Q 


(17) 


where 3X/3Q is a 7 by 11 matrix with initial conditions given by 


-12L s I, and M . o 
3X 0i 


(18) 


where I is a 7 by 7 unit matrix. For convenience, the following matrix 
notations are introduced: Let 


3X 

9Q 


m.A^.c 

_ M. - 

4 . 

M ] 

L v J ax 

3Q 

\ 3X 01 

3q / 


(19) 
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where 


U = — — = 7 by 7 state transition matrix 

dX n: 


V = = 7 by 4 control transition matrix 

9q 


A = 7 by 7 matrix 


C = 7 by 1 1 matrix 


1. Computation of the A matrix. From Eqs. (15), (1), and (2), 


A = 


3p7r T p 

~ I— T 


+ 4< i/m 
l/m 


0 j 3 


1 F T . 3 
m Ii / 


where the dimensions of submatrices are as indicated. (8£/9r) is obtained 


by differentiation of Eq. (13), and fromEqs. (11), (12), and (14). 
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2. Computation of the C matrix 
a. Derivation of (9#/ 9X Q .). From Eq. (15), 


d£_ 

a *oi 


ffn 

ar oi 


9H. 
1_ 

97 0i 


8F Ii 

a7 0i 


97 0i 


(21) 


where 


3F t . 

Ii 


3f. 

l 


3(r 0i , v 0 .) 9(r w , v„.) 


m 


From Eqs. (9), (7), (8), (5), and (4), 


9f. 


8(r 0i' v 0i> 


= f. 
1 


i an i ac 


n ap™ c 9p^_ 

op op 


2!°R + .I 

<>P. P. 


dP. 


8(r 0i’ 7 0i> 


where, from Eqs. (7) and (8), one can show that 


9p op _ _^i_ 8r l _ , . . 9c 

9p N ’ 9p (n i + 2l 2 p op )> 9 p 

a o op r op 


,C 1 + 2c 2 p op» 
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where 
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9p p 

- a 


dv. N 
1 o 


3H i _ 2,ip a J 1 3.1 2 3c 1 8p op 

8'i ' c 2 ‘ I n 3p op " c ‘PoplT 


From Eqs. (15), (2), (9), (8), (7), and (4) 


ajr 


0 


0 

a? Ii 


3L ? 

ap 0 


P 0 m 

9H. 
1 


9H. 
1_ 

. 8P ° . 


3p 0 

m 


(24) 


where 


!h . |_i_ , l an 8p °p i 

1 | P a n 8 Pop 3 Pa c 


9Pr 


9c 

P 


3 Ppp I 8p a 


op 


9p 


0 


8P. 

a Pr 


■ 


(i+3)/2 


i=l 


E. Search Algorithms 


1. Proble m statement . Let the dependent variables be 

Y(t F ) = X(t F ) - X B 


(25) 
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where = (x^, v^, 0) is a seven vector representing spacecraft position 

and velocity bias to the target planet. Let the target state vector be 
where Y = [x (t_), v (t„), mft^)] is the target planet position and velocity 

W p J? p Jc r 

vectors and the desired spacecraft final mass, if any, at final time, tp. 
Since Y(tp) is a function of independent search variables Q, where Q is a 
subset of all available independent variables, (r'(tQ), v(tg), m^), ffj, Pj, 
v. , Vg, tp] , in a linear approximation, the solution to the following 
equations gives the required corrections 6Q to the independent variables. 


-2X . 6Q = (Y - Y) (26) 

9Q w 


In actual nonlinear problems, procedure Eq. (26) is performed iteratively 
many times, until a satisfactory solution, (Y - Y ) = 0, is attained. 


2. Derivation of (9Y/9Q) . Due to the discontinuities arising in the 
thrust controls, the equations of motion and the variational equations are 
integrated piecewise with reinitialization performed at each discontinuity. 
Propagation and accumulation of the partial derivatives to the final time is 
required to obtain (9Y/9Q). 


a. Propagation of the state transition matrix. To obtain 
[8X(tp)/9X(L)j one must propagate stagewise information using the following 
chain rule: 


9X(tp) 9X(tp) 8X(tp_j) ax(t. + 1 ) 

1 — - ■ - ■ ■ ■■■■ • • • i i — i .i 

9X(t.) 9X(tp j) 9X(tp_ 2 ) 9X(t.) 


(27) 
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where 


ax(t.) 


ax(t. 


i+r 


ax 


Oi 


is obtained in the process of piecewise integration of the variational equations. 
Therefore, [aY/a(x Q , v Q , in^] is obtained in this manner with t. = t Q . 

b. Propagation of control matrix for control component v.. 

Since represents control applied only during t, to t i+1> [3X(t)/3v i ] * 0 
only for t. < t < t. + j. To propagate this control effect to the final state, the 
following computation is required: 


ax(t F ) ax(t F ) ax(t. +1 ) 


(28) 


3v. 

i 


3xit i+i> »“i 


where [8X(t. +1 )/8v. j is available at the end of the variational equation 
integration for the i-th power stage. 

c. Propagation and ac cumulation of control matrix for control 
components a and (3 r Since and Pj represent controls applied during 
tj to t I+1 , [aX(t)/a(a r Pj)] t 0 only for tj < t < t J+r 

Since the angle stages are designed to be larger than or equal to the 
power stages, each Atj contain several Ah' s. Therefore, accumulations 
and propagations of the piecewise control matrix must be performed. 
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Let tj = and tj. + ^ = implying Atj contains jf-power stages of equal 

duration. Then 


8X(t 

da 


i+r 


> ax(t I+1 ) ax(t k+n ) 

8X|t k + n) 8 “l 


( 29 ) 


Propagation of Eq. (29) to the final time as in Eq. (28) leads to 



ax(t F ) ax(t I+1 ) 

ax(t I+1 ) a«j 


The same procedure applies to obtain (3Y/3Pj). 


(30) 


d. Accumulation of control matrix for control component p^. Since 

Pq is a control applied all through the flight duration, procedure (29) is used 
— ► 

to obtain (8Y/9 Pq). Let k be the total number of power stages. Then 


e. 


8Y 

9p 0 


■£ 


Derivation of (9Y/9V-J. 

±5 


ax(t F ) 8X(t n ) 

From Eq. (25), 



(31) 


(32) 
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f. 


Derivation of (9Y/9tp). From Eq. (25), 


3Y £<V) aX B . , 8Y w (t F> 

®*p ^p ®*p 

3. Numerical search technique. In this section, a linear analysis 
technique used to solve for 6Q of Eq. (26) is given. Before proceeding to 
the detailed discussion of the analysis, row and column scalings will be 
performed to normalize Eq. (26). Let the dimensions of dependent variable 
space be n Y and that of independent variables be n^. Generally n^ > n^. 
Scaling factors for the independent variables (S^, i = 1, .... n^) are the 

conjectured largest step sizes within which the linear approximation 
[Y(Q + 6Q) = Y(Q) + (9Y/9Q) 6Q J holds. Scaling factors for the 

dependent variables (T\, i = 1 n Y ) are the accepted tolerance of the 

dependent variable deviations from the desired value. With these scalings, 
Eq. (26) can be transformed into a normalized form 

mx = y (34) 


where 


X. * 6 Q. / S . , i 1, 2, ***i 

ill W 

y. = (Y . - Y.)/ T. , j = 1, 2, • • • , n v 
wj y j J Y 

m. . = 9y./9x. 
ij Jl J 

Here, bold face letters are used exclusively to denote matrices, including 
raw and column vectors. In this normalized expression, convergence is 
considered to be attained if || y || ^ 1. The linear neighborhood constraint 
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on independent variables is implied by 


2 

x. . 
1 


1 X || < 1, where || X | = £ 

i 

Solution of Eq. (34) is obtained first by performing a singular 
value decomposition of m. The information gained in this singular value 
analysis is used further to control the selection of a particular solution. 
These policies of solution selection are of two types. The first type, called 
"minimal control policy," attempts to solve Eq. (34) using minimum || x || . 
The second type, called "final mass optimizing policy," attempts to solve 
Eq. (34) while maximizing the final mass. 

a. Singular value decomposition of matrix m. Let m be a real 
n Y by Oq matrix; there then exist matrices U, S, and V such that 

m = USV T (35) 

where u and V are square orthonormal matrices of orders n Y and nQ 

respectively (i.e. , UU^ = U^U = VV^ = v^V = I). 

/ 

S is a n Y by n^ matrix where only nonzero elements are on the 
principal diagonal (singular values) (i.e. , s.. = 0 for i j) and 


S 


11 


> s 


22 


> S 


33 


> S 


n Y n y 


(36) 


The basic algorithms to compute U, V, and S were given in Ref. 5 and the 
program that performs this computation exists in the JPL computer library. 


Consider the following orthogonal coordinate transformations of 
vector spaces X and y into x' and y': 

X' = V T x, y' = U T y (37) 
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Matrix m, in the representation of x' , y' coordinate systems, would be 
rn' = U T rnv o U T USV T V = S, and Eq. (34) is reduced to the following 
simple form: 


sx' = y' (38) 

Since S is nonzero only along the principal diagonal, i. e. , 



then Eq. (38), can be written as two separate equations 

SjX' = y' and = 0 
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The most general solution of Eq. (38) can be expressed by 



( 41 ) 


where 



Since s ii = 101 • c can be any arbitrary vector of dimension (n^ - n^), 

orthogonal to first m components of X' space. 


As one must be aware, n^ > implies that the problem posed is 
underdeterministic. The multitude of possible solutions given in Eq. (41) 
merely indicates this fact. 

b. Method of solution selection for minimal control policy. 

By definition, minimal control policy implies that Eq. (34) is satisfied while 
minimizing ||x|| (call it minimal length). Note that ||x|| = ||x' || holds due to 
the orthonormality of v used in Eq. (37). 
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According to Eq. (41), | X* | = || XJ. || + || X ^ || , since X'j is 

constrained by Eq. (41), minimum length J| X | is obtained if one set Xjj = C 
= 0. Therefore, the basic choice of solution for Eq. (38) is: 



n 


Y 


X 



n 


Q 


If one performs inverse transformation of Eq. (37), 



Sj'V 


v s 

x = vx' = V 

0 

= V 

0 




. 


(42) 


i. e. , 



(43) 
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Let Z - 3m(tf)/3x be a 1 by row vector, and m(t^) is the 
spacecraft final mass. In the x' representation Z is transformed to 


9m(t,) 

Z = 7 =- = ZV (45) 

8X 


Then mft^.) = mg(t^) + z'x' , where m^t^) is the current mass, and m(tj) is 
the linear estimate of the mass after corrections x' are applied to the 
independent variables. Further, let 


Z' 


n Y n Q ” n y 


Then 


m(t f ) = m Q (t f ) + zjx| + zJjX^ 


(46) 


As was shown in Eq. (41), the addition of an arbitrary vector x^, which is 
orthogonal to Xj, does not disturb the constraint s'x' = y 1 or equivalently 
mx = y. Therefore, it is possible to construct a vector that will modify 
mft^) by making Z^jXjj as large as possible. Since ZjjXjj is linear and 
unbounded, a notion of maximum is not valid. However, since || x ijjj < 1 is 
also necessary in this type of linear iteration procedure, we can assert that 
for ||Xjj|| = 1 , the maximum expected z jj x jj can be obtained if one chooses 
X II </ { z \j) • w itt l length 1; i.e. , 
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The basic solution given in Eq. (43) is generally not directly applicable 


for highly nonlinear problems. One often would encounter the cases where 
even the minimal length solution of Eq. (42) has length || X || exceeding 1. 
This violates the nominal linear domain constraint and is undesirable. To 
resolve this dilemma, a careful inspection of singular values Sjj, 


S 


. . . , s is in order. Since X' = S *y', X’ can become large if 

L.t. i 


some components of 8^ becomes rather small. As are given in the order 

of decreasing magnitude, one can readily examine such situations. When the 

ratios of Sj j to (defined as condition number in Ref. 4) becomes larger 

5 

than some number (e. g. , ~ 10 ) for i k + 1, one may consider that the 

given matrix m of Eq. (35) actually is ill-conditioned. If this situation is 
encountered, it is likely that one is dealing with a nearly correlated matrix 
m whose effective rank is k (< ny) instead of ny. Then, the proposed 
solution is obtained by replacing Eq. (43) by 


(44) 


In the event that the condition number of matrix m is not large, uniform 
scaling of the solution given by Eq. (43) is recommended to restrict 

INI < i- 



c. Method of solution selection for the final mass optimization. 
If one were concerned with the outcome of the final mass, and wished to 
maximize the mass while satisfying the constraint mx = y, the following 
scheme is suggested. 
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( 47 ) 



Then 


Z II X II “ 


V 1 Z II 



(48) 


This is the mass gain one can expect if one modifies the minimal length 
solution by the addition of component Xjj. 

Therefore the best mass optimizing solution is 


1 . e. , 


X = VX = V 







— 

x i 


s l" 

‘y’ 


s^uy 


= V 

— 


= V 


X' 


(Z* 

n ,T 


< ZT »n 

II 

» _ 


vi 

l z 'nll 


VI^SII 
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— » 


n 


Y 


X i ^j V ij S.. + ( 

j=i JJ n=i \ 
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Note here that the index can take values k < n^, when the 
matrix given is ill-conditioned. 


d. Comments on the search techniques. In principle, if the search 
variables are updated iteratively using Eq. (50) for corrections, it 
eventually will satisfy the boundary conditions where || y|| < 1 is reached. 
Further, if the mass increment indicator VIM of Eq. (48) becomes 
smaller than the pre-assigned number, e. g. , 1 to 2 kg, one may consider 
that an optimal final mass is attained. To date, extensive tests using 
minimal control policy solutions given in Subsection E-3-b have been per- 
formed. The results are very satisfactory in most cases. The algorithm 
described in Subsection E-3-c is still under investigation, however it 
appears promising. 


V. FUTURE DEVELOPMENT 

It is intended that future work will include: 

(1) A thorough test of the optimizing algorithm. 

(2) Further refinements and verification of the representations of 
low-thrust subsystem characteristics. 

(3) A more detailed development of the requirements in the 
interface with the OD and guidance programs, particularly the 
input/ output specifications. 

(4) Investigations into the modeling of the lower order perturbations 
and verification of current findings that these perturbations will 
not be a major problem. 
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The final product from these studies will be a set of software requirements 

specifications for a flight quality trajectory program. 
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